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VARIANCE REDUCTION IN MONTE CARLO ANALYSIS OF RAREFIED GAS DIFFUSION 

by Morris Perlmutter 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 

The present analysis uses the Monte Carlo method to solve the 
problem of rarefied diffusion between parallel walls. The diffusing 
molecules are evaporated or emitted from one of two parallel walls and 
diffuse through another molecular species (see fig. 1). The M.C. anal- 
ysis treats the diffusing molecule as undergoing a Markov random walk 
and the local macroscopic properties are found as the expected value of 
the random variable, the random walk payoff. By biasing the transition 
probabilities and changing the collision payoffs we can retain the ex- 
pected Markov walk payoff but reduce its variance so that the M.C. re- 
sult will have a much smaller error. 

As shown in figure 1 the Markov random walk of the diffusing 
molecule can be represented by the sequence (X 0 ,X 1 ,X 2 , . . . } . X n refers 
to a point in velocity and position space (-V n ,Z n ) taken by the diffusing 
molecule immediately after the n th collision. The probability density 
corresponding to this random walfc is given by 

H(X 0 ,X 1 ,...)dX 0 <%... = E 0 (X 0 )K(X 1 |X 0 )K(X 2 |X 1 ),...dX 0 <%.. . (l) 

The birth distribution E Q (X 0 ) refers to the probability of the 
molecule originating at X Q , Assuming the molecules leaving the wall 
are in thermal equilibrium the birth distribution for the z component 
of the dimensionless velocity of the molecules leaving the wall is 
given by E Q ( v Zo’ Z ~ 0) = 2v Zo exp(-v| 0 ). The transition probability 
K(Xn+i| X n )dX n+ i can be written as the product of a transport proba- 
bility T(Z n -» Z n+1 | v n )dZ n+1 and a collision probability 

C(v n -* v n +l| z n+l) dv n+l- The transport kernel T gives the probability 
of leaving Z n and reaching Z n+1 at the n + 1 collision while the 
collision kernel C gives the probability of a molecule at velocity 
v n reaching a new velocity v n+1 after the n + 1 collision. For the 
present model the transport kernel can be written in dimensionless form 
as T(Z n - Z n+1 |v n ) =(V|v Zn |)exp(-|Z n+1 - Z n |/|v Zn |). The collision 

kernel, assuming the molecules come out of collision with a Maxwellian 
distribution, will then not be a function of the previous velocity or 
position and can be written as C(v Zn ) = 0-/ aJtt) exp (-v| ) . We can write 
the event payoff after the nth collision. as P(X n ). Then the payoff 
for the random walk t) is given by 

O oo 

V ' E P(X n>- 

iteO 

The expected value of the random walk payoff is given by 
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A = <T1 0 > 




P(Xn) 


E 0 (X 0 )K(X 1 |X 0 )K(X 2 |X 1 )...dX 0 


• (2) 


The payoff term after each event can he written as . P(X n ) = 
P( v n) T ( Z n "* z sl v n)' The P( y ) is some func,fcion of velocity of the 
molecuie after- a collision and t(Z -» Z g |v) is the scoring probability, 
the probability of the sample molecule reaching the scoring position 
Z s from the last collision position without having another molecular 
collision. We can thus write (t] 0 ) = n s (p(v)v z ) s /p 0+ ; where n g is 
the local molecular density at the scoring cross section and p 0+ is 
the molecular flux leaving the emitting surface. If the macroscopic 
quantity desired at the scoring position Z g is the mass flux, |j. g , 
then p ( v) is given by p^ = ±1 for v^ ^ 0. Similarly, if the local 
molecular number density n g is desired at Z g we have p n = ± l/v z 
for v z ^ 0. The scoring probability can be found to be r(Z n -» Zg|v Zn ) 
= exp(-|Z n - Z s |/|v Zn |). 


The potential payoff of a molecule leaving X i is given by 
T] i = P( x i) + P( x i+i) + ••• • Averaging this we can write the expected 
potential payoff conditional to leaving Xj_ as 


/ ^ 


»l(Xl> * <1i|Xi> = / ^ 

J t=0 


P(X 1H )K t (X i+i |Xi)dXi H 


= P(X t ) + J W i+1 (X i+1 )K(X i+1 |X i )dX i+1 (3) 


where 


x ^,( x i+^| x i) - J* • • • J' ^ x i+l|A-i) • • • x ( x i+^| Xi+^-iJ^i+i" . *dX^ + ^_2_. 

We can similarly calculate the expected value of the square of the 
potential payoff of a molecule leaving X.^ as 


Q i (X i )= (nS|Xi> = 2P(X i )W i (X i ) - P 2 (X i ) + j r Q i+1 (X i+1 )K(X i+1 |X i )dX i+1 


if 

n^O J 


f 2P ( X i + n>'W ) W -?& 


x+n 


(4) 


The expected variance can now be written as . 

i|i 0 >wy*l- <%> 2 r;/«o(Xo) E o(Xo) d Xo - <io> 2 - ' 

In the M. C. process we pick a birth velocity randomly from E 0 (v Q ), 
then score P(X Q ). We then randomly choose the position of first colli- 
sion from K(X 1 |X 0 ). We then score P(X 1 ). The particle history is 
continued until it is incident on one of the walls. The payoff for the 
random walk is then given by t) q . This process is repeated for W sam- 
ples. -The expected payoff is then given by 
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whose value is the desired macroscopic result appropriately representing 
the molecular flow rate or molecular density. The 95 percent confidence 
interval of rj is given by 

f(i 0 ) 

|(, o) < 6 = 1.9S-^- 

where cr(q 0 ) is the standard deviation of r) 0 and can be obtained from 
the M.C. calculation by 

° 2 (ti o) ** ( N - i ~ Y,' ^oi ' (^o ) 2 
i=l 


The analog M.C. calculation is the case in which the sample payoffs 
are given as p(v) and are scored only when passing the scoring position. 
In table I the molecular density of the diffusing molecules at the scoring 
position Z„ = l normalized by the density of the molecular flux emitted 
at Z = 0 is given under the heading "Analog." Also given are the re- 
sults where the payoff is p(v n )T(Z n -» Z s |'v n ) = P(v n ,Z n ). This is scored 
after each collision. These results are given under the heading "next 
event . " 


We wish to reduce the variance , a (q 0 )'. sob that the- error in. the M.C. ■ 
sampling will be smaller.' To do this we bias the probabilities as follows 
The probability of the- biased walk is now given by H*(X 0 ,X 1? . , . ) = 


E*(X )K*(X 1 |X )K*(X 2 |X 1 ). . . . If we then distort the payoff to be 

„ p/v , vioi + P(x , vviaiv + 

1 ~ ° V x o> 1 EJ(X 0 )K*(XjJX 0 ) 


We can see that the expected value of the random walk payoff is given by 
dX Q dX-|_... = X and is unchanged by the biasing process. 

A useful simplification is to write K*( x n +ll X n) = ^(^n+ll x n) {-'-n+l^n+l V 
I (Xn)]. The biasing function must satisfy the usual probability normal- 
ization requirements. The probability density corresponding to the biased 
random walk is given by 

r ii( x i)l 

[E 0 (X 0 )I 0 (X 0 )]^K(X 1 |X 0 ) 


H = 


Wj 


k(x 2 |x 1 ) 


i 2 (x 2 ) 




while the payoff for the random walk becomes 

P(_X 0 ) P(X X ) P(X 2 ) 

'o 


< = 


I (X ) 
o v o' 


x l( x l) 


i,(x,) 


For the biased case, the expected potential payoff for the random walk 
is given by W*(Xj_) = [W i (X i )/l i (X i ) ] . 

Similarly, we can find the expected value of the potential square 
payoff for the biased random walk as 



4 


K n*' X i+nl X i^ 


y ^i+n 


00 y-» 

i( x i) = x . (X- ) X/ / ^^ X i+n^i+n( X i+n) “ ^( X i+n^ ^i+n^ X i 

1 1 n=0 J 1 " !l 1 ' “ 

Since (tj* 2 ) = J' E*(X o )Q*(Xo)dX o we can write the biased - expected 


payoff as 


r C e (x ) 

<< 2> = J iy^T [ 2P ( X 0> M 0< X C,> -^(Xo) ]<3X 0 + J ^lyBPtxpw^xp 

- ^(xpldXL * j Y^r y prCx^^) - r 2 (x 2 ))dx 2 + ... 


( 5 ) 


We wish to minimize (ti* 2 ) with respect to I n ( x n ) subject to the condi- 
tion that 

I. 


y ’... y ' t E 0 ( X 0 ) l 0 ( X 0 ),[ K ( X 1 | X 1 , ^ y ] [ k ( X2 | Xi ) ^ y ] ... 

-/ E n( x n) I n( x n) ax n = 1 
Using the calculus of variations we find 

[2P(X n )W n (X n ) - F 2 (X n )] 1 / 2 


dX„ dXn 


I n (Xn) = 


( 6 ) 


( 7 ) 


y , E n (X n )[2P(X n )W n (X n ) - P 2 ( X n )] 1 / 2 dJ^ 

Since 2P(X n )W n (X n ) - P 2 (X n ) = Q n ( x n ) “ J^'^n+l^ X n+l^ K ^ X n+ll X n^ ( ^ X n’ this 
result implies that the importance function for the kernel biasing is 
proportional on the contribution to the variance of the molecules coming 
out of collision at X n but not the future contribution to the variance 
from the later collisions. If we neglect higher order collision terms 
we can approximate the bias function by 

-w - P(Xn) 


/VVPOgaXn 


(a) 


Now the transition kernels are biased so that more samples are taken 
with values that give larger contributions to the payoff. The biased 
payoff is now given by 

p(xj r 

P*(XJ = y-2-y = yE n (X n )P(X n )dx n = A n (9) 


The payoffs are constant values for each .event. This will reduce the 
variance of the random walk payoff. 

The birth payoff is given by P*(X Q ) = A Q = J* E 0 (X 0 )P(X 0 )dX Q . This 
integral can be evaluated by numerical integration. Then a birth velocity 
can be randomly picked as usual and the unbiased procedure continued as 
before to find the next payoff P(X-j_). The result for this procedure is 
shown in table I labeled birth bias for the same cases as before. 
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We can continue the biasing to. the first collision, P*(X 1 ) = = 

yE 0 (X 0 )K(X 1 |X 0 )P(X 1 )dX 0 dx 1 . However, it is difficult to evaluate 
by numerical integration. Instead we can approximate the importance 
function by 

- P(Xi) 

- T 0 (Z 0 ,v 0 ) CP( Zl ) 

where T 0 (Xq) is given by J T(Z q -*• Z 1 |v Q )dZ 1 and CP(Z 1 ) = 

“ C(v 1 )P(v 1 ,Z 1 )dv 1 . The payoff is then given by P*(X-]_) = 

T 0 (V Z o)^ Z l)* To evaluate ' this we randomly choose v Q from E Q (v o ) 
and Z 1 from T(Z Q -> Z 1 | v o )/T 0 (X Q ) . With values for v Q and Z ± we 
can then evaluate the first collision payoff P*(X^). We can then con- 
tinue the unbiased procedure retaining the v q and Z^. However since 
the Z± was found from the biased distribution the future payoffs must 
be weighted by P*(X n ) = P(X n )T 0 (X 0 ). These results are shown in 
table I under +1 term bias. 

This process can be continued to the second collision using 

i 2 (Xg) = p(x 2 )[t o (x o )t 1 (x 1 )cp(z 2 )]- 1 (11) 

The T-^X-l) is given by [ Z]j (Z ± - Zg | v 1 )dZ 2 where 2^ is L or 0 

^1 

depending on v 1 ^ 0. The payoff is then given by P*(Xg) = 
T(X 0 )T(X 1 )CP(Z 2 ). To evaluate P*(Xg) we already have v Q and Z ± , we 
evaluate v g from C(v 2 ) and Zg from T^ -* Zg |v 1 )/T 1 (X- l ) . We can 

then evaluate the P*(Xg) payoff and we can continue either in the 
biased or unbiased fashion in this manner. The results for the biasing 
for three terms is shown in table I under +3 term bias. 

Finally we can bias continuously, however in this case the sample 
history will not end because of the transport kernel biasing. The 
sample histories were then ended using Russian Roulette, samples were 
followed until the weighting T(X q )t(X 1 ) . . . was less than 0.001, then 
if a rando ml y picked number, uniform between 0 and l,was greater than 0.1 
the history was ended, if less the sample weight was multiplied by a 
factor of 10 and the process continued. These results are in table I 
under Russian Roulette. The results indicate significant reductions 
in variance can be obtained using biased sampling techniques. 



TAELE I. - NUMERICAL RESULTS 


Sampling 

Density 

Mean 

Time, 

Sampling 

Density 

Mean 


ratio, 

deviation, 

min 


ratio, 

deviation, 


n / n o+ 

a 




0 


10 

000 Samples 


10 

,000 Samples 


L/A = 0.1 



L./A = 10 

Analog 

0.848 

0.923 

.0.13 

Analog 

0.1202 

0.578 

Next event 

.844 

.867 

.15 

Next event 

.11965 

.527 

Birth bias 

.845 

.741 

.12 

Birth bias 

.11965 

.527 

+1 Term bias 

.850 

.1568 

.42 

+1 Term bias 

.1171 

. 463 

+3 Term bias 

.855 

.116 

.64 

+3 Term bias 

.1195 

.402 

Russian Roulette 

.854 . 

.115 

.61 

Russian Roulette 

.1172 

.244 


L/A = 1 



L/A = 50 

Analog 

0 . 4914 

— 

0.7494 

0.31 

Analog 

0.0261 

0.224 

Next event 

.4881 

.7188 

.35 

Next event 

.0263 

' .238 

Birth bias 

.4881 

.7185 

.32 

Birth bias 

.0263 

.238 

+1 Term bias 

.4958 

.5728 

.59 

+1 Term bias 

.026.3 

. .2759 

+3 Term bias 

.487 

.336 

.95 

+3 Term bias 

.0246 

.183 

Russian Roulette 

.491 

.259 

2.33 

Russian Roulette 

.0245 

.168 


Time. 

min 


2.14 

2.33 

2.31 

2.48 

3.22 

13.97 


10 . 78 
11.61 
11.60 
11.34 
14.17 
39.50 


<<r 
O' 

• NONDIFFUSING SPECIES 

I O DIFFUSING SPECIES 

tol SCORING CROSS SECTIONS, Z 



Figure 1. - Analytical model. 



















